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Simulation study of pressure and temperature dependence of the negative thermal 

expansion in Zn(CN) 2 
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Pressure and temperature dependence of the negative thermal expansion in Zn(CN)2 is fully in- 
vestigated using molecular dynamics simulations with a built potential model. The advantage of this 
study allows us to reproduce all the exotic behaviours of the material, including the negative thermal 
expansion (NTE), the reduction of NTE with elevated temperature, the pressure enhancement of 
NTE and the pressure-induced softening. Results of the study provide us detailed data to link the 
properties in the energy space and the real space, giving us insights to understand the properties 
and the connections between them. 



PACS numbers: 65.40.-b, 65.40.De, 63.20.Dj, 02.40.-k 



I. INTRODUCTION 



Negative thermal expansion (NTE) is a rare and 
counter-intuitive phenomenon found primarily in low 
density materials with crystal structures that are net- 
works of linked coordination polyhedra. The study of 
these materials is not only of fundamental scientific im- 
portance, but also has many technological applications 
such as aerospace technologies^-, optica^, and electronics^]. 
While much effort has been put into finding new mate- 
rials and investigating the origin of NTE in them, much 
less attention has been paid to the change in NTE be- 
haviour subject to heating and stress, which holds great 
importance for the possible applications of the material. 
For example, due to stresses and heating, problems such 
as phase transitions of the NTE filler and thermal expan- 
sion misfit between the NTE filler and matrix are always 
encountered i n d esigned composites with tailored ther- 
mal expansiorpEl. 

In this paper, we conduct a simulation study of 
Zn(CN)2 focusing on the pressure and temperature ef- 
fects on its negative thermal expansion. We chose this 
material for several reasons. Firstly, Zn(CN) 2 is a well- 
known representative NTE materiaP. It has a frame- 
work structure consisting of tetrahedral groups of atoms 
linked by diatomic rods of C-N and has exceptionally 
large isotropic NTE of a lincar = —16.9 MK _1 (twice as 
large as that of ZrW20fl). Secondly, the material shows 
a variety of exotic properties in experiments^, including 
reduction of its NTE on heating, pressure-enhanced ther- 
mal contraction, and pressure-induced softening, none of 
which are fully understood. Thirdly, although previous 
DFT calculations of Zn(CN) JaHHUI can explain the origin 
of NTE of the material in terms of Griineisen theory, the 
calculations do not draw a clear link between the val- 
ues of the Griineisen parameters in energy space and the 
structural vibrations in real space. Moreover, the calcula- 
tions have not been able to reproduce the aforementioned 
exotic properties of the material observed in the experi- 
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FIG. 1: Difference plot of the electrostatic potential arising 
from the rank 4 DMA and rank MulfhP^] charge models vi- 
sualized on the van der Waals surface of (a) a Zn-C-N ordered 
cluster, (b) a Zn-N-C ordered cluster. The colour bar shows 
values of the electrostatic potential difference in kj/mol. The 
average percentage error in the electrostatic potential made 
by the point charge model on this surface is 1%. 



ments. 

Here, we have built a Zn(CN) 2 potential model based 
on first-principles calculations. Lattice-dynamic calcula- 
tions and large-scale molecular dynamics (MD) simula- 
tions were carried out for the material using this model 
which was justified by comparing against the available 
experimental data. The results in both energy and real 
space provide us fundamental clues to understand the 
NTE as well as the related exotic behaviours of Zn(CN) 2 . 



II. BUILDING THE MODEL 

We started with total-energy calculations for a 
[Zn(CN)4] 2_ cluster. The —2 charge comes from the fact 
that every zinc is shared by four neighbour atoms (C 
or N), and that the bond should be highly ionic accord- 
ing to an initial judgment of zinc having much smaller 
electronegativity than carbon and nitrogen. Four CN _ 
molecular ligands around each Zn 2+ would give us a —2 
charge on the cluster. Geometry optimization of the clus- 



TABLE I: Potentials used for various interactions in the model. E is the energy, r is the distance between the two atoms, ro 
is the equilibrium distance, 9 is the bond bending angle of C/N-Zn-N/C in degrees, 9q is its equilibrium value and <p is the 
linear-three-body angle of Zn-C/N-N/C. Units are eV for energy and A for length. The prefactor of the harmonic-bond-bending 
term has unit of eV/rad. 



Potential 



Form of the Potential 



Type of Bond 



Values of Parameters 



Morse potential 



E M = D{[1 - exp (-o(r - r„))] 2 - 1} 



Harmonic three-body potential 



Linear three-body potential 



Buckingham potential 



E H = (1/2)K{9 - Oof 
E L = K(l-coB<p) 
E B = Aexp(-r/p)-(C/r 



Zn-C 


D = 


= 0.2432, a = 2.917, r = 2.109 


Zn-N 


D = 


= 0.1795, a = 2.701, r = 2.157 


C-Zn-C 
C-Zn-N 

N-Zn-N 




K = 1.5157, 6 = 109.47 
K = 1.2695, 6>o = 109.47 
K = 1.0233, = 109.47 


Zn-C-N 




K = 1.1968 


Zn-N-C 




K = 0.6359 


C-C 

C-N 

N-N 


A = 
A = 
A = 


= 2806.9, p = 0.2667, C = 17.67 
= 2365.0, p = 0.2825, C = 20.88 
= 1992.7, p = 0.2874, C = 24.67 



ter resulted in a perfect tetrahedral conformation. 

Total energies for different configurations of the clus- 
ter (with bond stretching and angular distortions) were 
compute d usin g DFT in GAMESS(US) with the PBE0 
functiona l 14 * 15 ^. Correlation-consistent basis sets up to 
aug-cc-pVQZ were tested, and an aug-cc-pVTZ basis set 
was found to have sufficient accuracy and no significant 
basis set superposition errors (BSSE) at different config- 
urations. 

Various interatomic potential forms were then fit- 
ted to the calculated energy curves to obtain the ini- 
tial potential parameters. The short-range interactions 
are the Morse potential for Zn-C/N, a harmonic three- 
body-bond-bending term for C/N-Zn-N/C and a linear- 
three-body term for the angular distortion of Zn-C/N- 
N/C. The long-range van der Waals interactions are de- 
scribed by a Buckingham potential with parameters from 
WmiamsP. 

Two types of clusters with Zn-C-N order and Zn-N-C 
order are used. The multipoles on each cluster were cal- 
culated by distributed multipole analysis (DMA)P3 using 
CamCASP 1 -^. The effective point charges on the atoms 
were then obtained by fitting to the electrostatic poten- 
tial from the rank 4 (hexadecapo le) d istributed multi- 
poles using the MULFIT prograrri 19 l 20 l Figure [I] shows 
the difference in the electrostatic potential between the 
point-charge model and the DMA result. The root mean 
square of the difference is less than 4 kJ/mol and 8 
kJ/mol for clusters of Zn-C-N and Zn-N-C, respectively, 
corresponding to about 1% relative difference in the elec- 
trostatic potential around the clusters. The averaged ef- 
fective point charges on each atom are (in electron unit) 
+1.14 for Zn, -0.21 for C and -0.36 for N. 

The initial potential model was then refined by refit- 
ting to DFT energy surfaces of various cluster configura- 
tions. The point charges were fixed during this process. 

This gave us final potential parameters that implicitly 
incorporate the effects of higher-ranking multipole mo- 
ments and atomic polarization. Due to the high strength 



of the C=N bondMEU, ^ mg g r0U p was treated as a rigid 
rod in all cases. The potentials with their parameters are 
listed in Table |U 



III. COMPUTATIONAL METHODS 

Harmonic lattice dynamics (HLD) and quasi-harmonic 
lattice dynamics (QHLD) calculations were carried out 
using GULF^l 

Molecular dynamics (MD) simulations were performed 
using DL_POL\^for a 10 x 10 x 10 supercell containing 
10000 atoms with periodic-boundary conditions. A con- 
stant stress constant temperature (NaT) ensemble with 
a Nose-Hoover thermostat 2 - was used. The long-range 
Coulomb interactions were calculated using the Ewald 
method with precision of 10~ 6 . The equations of motion 
were integrated using the leapfrog algorithm with a time 
step of 0.001 ps. A total of 20000 time steps were used to 
achieve equilibration. At different temperatures and pres- 
sures, snapshots of atomic trajectories after equilibration 
were recorded every 0.02 ps up to a total of 50 ps for the 
follow-up analysis. A disordered model of Zn(CN) 2 was 
constructed by randomly switching C and N atoms in the 
supercell. 



IV. AN INITIAL TEST OF THE MODEL 

We compared some basic quantities obtained from the 
model against experiment. The optimized structure gave 
the cell parameter as a = 5.9176 A, and the Zn-C/N 
bond length as r = 1.9782 A compared to a = 5.9227(1) 
A and r = 1.9697(3) A from an experiment at 14 K 24 . 
The model is found to be stable in lattice-dynamic cal- 
culations, and the calculated phonon frequencies are in 
good agreement with the spectroscopy dataP^, as shown 
in Table [H] The dispersion curves shown in black in 




R M 

Wave vector 



R 



FIG. 2: Calculated dispersion curves of Zn(CN)2 along several 
high-symmetry directions in the BZ. Black is for zero pressure. 
The acoustic modes around M (0.5, 0.5, 0.0) and the middle 
point (0.25, 0.25, 0.25) of F-R are the first to become unstable 
beyond 1.0 GPa, as shown in light red. 



Figure [2] are in good agreement with previous DFT 
calculation j 10 1 13 1 . 



It is worth noting the close similarity between the dis- 
persion relations of wave vectors along T-X-R and R-M- 
r, as well as the near-mirror symmetry of the dispersion 
curves along T-R. This makes sense given the existence of 
the two interpenetrating cristoablite-like networks in the 
material. For the lowest-frequency mode at wave vector 
R, the two networks would translate like acoustic modes 
but out of phase with each other, leading to a dispersion 
relation that has the appearance of an acoustic mode but 
with a non-zero band gap. 
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FIG. 3: The calculated isotherms of Zn(CN) 2 from 0.0 GPa 
to 0.7 GPa with an 0.1 GPa increment (solid lines). The inset 
shows the calculated pressure enhancement in a compared 
to the experiment. The negative thermal expansion of the 
material will reduce on heating. 
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FIG. 4: NTE of the disordered system (in symbols) com- 
pared to that of the ordered system (in solid lines) indicating 
a trivial difference of the two systems. 



V. THERMODYNAMIC PROPERTIES 



TABLE II: Calculated phonon frequencies compared to the 
spectroscopy datcPM 'OM' stands for the ordered-model cal- 
culation. 'VC stands for virtual crystal approximation. 'Calc' 
and 'Exp.' represent calculated and experimental data, re- 
spectively. 'IR' and 'R' represents Infra-red and Raman ac- 
tive, respectively. 



OM Calc. 


173 


186 


328 


337 


463 


VC Calc. 


210 


220 


314 


349 


490 


Exp. 


178 


216 


334 


339 


461 


Irrep. 


Ti„ 


T 29 


E 9 


T 2s 


Ti u 


Type 


IR 


R 


R 


R 


IR 



A. Pressure and temperature dependence of the 

NTE 



The calculated NTE curves of Zn(CN) 2 under different 
pressures from 0.0 to 0.7 GPa with increments of 0.1 GPa 
are shown in Figure p3j The pressure-enhanced a (aver- 
aged over 50-300 K) at 0.0, 0.2 and 0.4 GPa are -12.62, 

— 14.09 and — 15.97MK -1 , respectively, compared to 
the experimental valuesPof -17.40(18), -18.39(27) and 

— 19.42(23) MK -1 at the corresponding pressures. The 
MD successfully captured the gradual reduction of NTE 
on heating which has been observed in X-ray scattering 8 . 
According to the MD, a v is -30.3 MK _1 at 300 K, much 
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FIG. 5: The pressure and temperature dependence of the 
bulk modulus of Zn(CN)2 calculated from the MD p-V data 
using thermodynamic relations. 



FIG. 6: The pressure dependence of the volume of Zn(CN)2 at 
different temperatures. The inset shows the change of volume 
with pressure up to 3.0 GPa at 300 K. A high-pressure phase 
settles down after the clear discontinuity at 2.1 GPa. 



lower than the value at 25 K, -47.3 MK" 1 . 

Experimental has confirmed that Zn(CN) 2 exits in a 
disordered form, i.e. the carbon and nitrogen atoms are 
randomly placed on their symmetric sites in the struc- 
ture. NTE of the disordered model from the MD calcu- 
ation is displayed as symbol plots in Figure [4] compared 
to that of the ordered model in solid lines. There is no 
significant difference between the NTE property of these 
two systems. 



B. Pressure and temperature dependence of the 
mechanical properties 

The bulk modulus is directly related to the coefficient 
of thermal expansion a v according to the classical form 



3i?7 

Yb 



(i) 



where R is the gas constant, the volume V and the bulk 
modulus B. The overall Gruneisen parameter is calcu- 
lated by summing over all modes with frequencies Wk,A 
and the volume strain s 



N 



1 Olu 



k.A 



k.A ' 



ds 



(2) 



At 300 K, the values of the bulk modulus Bq and its 
first derivative with respect to pressure B' Q at zero pres- 
sure obtained by using the 3rd-order Birch-Murnaghan 
equation of states to fit to the isotherms calculated from 
our MD are B = 34.47(31) GPa and B' = -4.2(5). Both 
are in good agreement with the experimental values of 
34.19(21) GPa and -6.0(7), respectively^. However, the 



values of the bul k mo dulus obtained from the previous 
DFT calculations^^ are much higher the experiment. 

We obtained the bulk modulus of Zn(CN) 2 at differ- 
ent temperatures and pressures by numerically comput- 
ing derivatives of the p-V data from the MD simulations. 
As shown in Figure [5j the bulk modulus doesn't change 
much with pressure, but changes enormously with tem- 
perature, namely it increases by about 40% on cooling 
from 300 K to 10 K. Referring to Figure |3j we found 
that a 5% decrease in volume caused by heating at zero 
pressure corresponds to as much as 60% decrease in the 
bulk modulus, while the same amount of volume decrease 
caused by compression would only reduce the bulk modu- 
lus by less than 5%. This means that the bulk modulus of 
the material not only depends on the volume change per 
se, but also on the means of changing the volume — by 
heating or compression. This breaks Birch's law of cor- 
responding states^lEIl The same anomaly has been ob- 
served experimentally in Zrx^Os, where the bulk mod- 
ulus increases by 40% on cooling from 300 K to K 33 . 

We found that the value of B' from a lattice-dynamic 
calculation using GULF 26 is 7.2 compared to its nega- 
tive value at 300 K. This suggests that all the mechanical 
contributions at T = are from the Zn-C/N bonds that 
become stiffened on compression. With elevated temper- 
ature, one can imagine that the Zn-C/N-N/C angle flex- 
ing starts to contribute to the change of volume on pres- 
sure. This mechanism costs much less energy than com- 
pressing the Zn-C/N bonds as suggested by values of the 
parameters in the Morse potential and the linear-three- 
body potential of the model (Table |l|. As a result, the 
bulk modulus decreases on compression hence the nega- 
tive B' , i.e. pressure- induced softening of the material. 
At high temperature, B' is expected to become less neg- 
ative on heating due to the rising energy cost of further 




• Experiment 

Calculated at 297 K 



\.S+«s., 



2 4 6 8 10 12 14 16 
Energy (THz) 




4 6 8 10 12 14 16 

Energy (THz) 



FIG. 7: Calculated phonon spectrum at 297 K compares to 
the experimerrtP^l. 



increasing the Zn-C/N-N/C angle vibrational amplitude 
on compression. Values of Bq at different temperatures 
are given in Table \UJ\ 



FIG. 8: The Phonon DoS calculated from the velocity auto- 
correlation function, (a) The overall DoS at 0.0 GPa. (b) The 
angular DoS of the C-N rods at 0.0 GPa. Arrows show that 
the change of peaks on heating. 



VI. PHASE TRANSITIONS 

In Figure [61 simulations were conducted for Zn(CN)2 
at different temperatures in a much broader range of 
pressure from —2.0 to 3.0 GPa, and there is clearly a 
phase transition of the material caused by compression at 
each temperature. With elevated temperature, the phase- 
transition pressure indicated by the discontinuity in the 
volume change increases, implying that the phase transi- 
tion may be triggered by some soft modes that can be sta- 
blized on heating due to the anharmonic term in their fre- 
quencies. Indeed, as shown in Figure [51 if we compare the 
dispersion curves calculated at a pressure beyond 1.0 GPa 
(in light red) to that at zero pressure (in black), we find 
softening of the acoustic modes around the zone bound- 
aries, especially the modes at M and the mid-point of 
T-R (0.25,0.25,0.25) which are the first ones to become 
unstable. The concurrent softening of the optic modes 
(~ 1.5 THz) directly above these acoustic modes suggests 
a possible hybridization between the acoustic modes and 
the optic modes, resulting in a k 2 energy behaviour^ at 
k / 0. 

The inset of Figure [6] shows the change of volume 
with pressure up to 3.0 GPa at 300 K. The volume dis- 
continuities at 1.2 and 2.1 GPa may suggest hysteresis 
in the phase transition. We found that the new high- 
pressure phase is orthorhombic with P2i2i2i space group 
(a = 10.72 A, b = 10.78 A, c = 10.88 A). This is 
compared to an orthorhombic phase with Pmc2\ space 
group found beyond 1.3 GPa in the X-ray diffraction 
experiment!^. 
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FIG. 9: (a) Total DoS at 300 K under different pressures, 
(b) Angular DoS of C-N rigid rod at 300 K under different 
pressures. Arrows show the change of peaks on compression. 



VII. PHONONS 

A. Density of states 

The Fourier transformation of the atomic velocity 
auto-correlation function (VACF) gives us the phonon 
density of states (DoS)^ of the material. To obtain the 
DoS at different pressures and temperatures, we used the 
trajectory data of atoms from MD to calculate the VACF 
of Zn(CN) 2 . 

The correlation function C(t) at time t = nAt can be 
expressed as, 
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FIG. 10: (a) Dispersion curves coloured according to the corresponding value of the mode Griineisen parameter 7k, a, where 
red corresponds negative values with colour strength representing values between —18 and (white), and blue representing 
positive values, 0-4. (b) Dispersion curves coloured according to projection of the eigenvectors onto RUMs, with red and blue 
contributions representing the extent of rotations and translations respectively, with colour strength representing values between 
100% and 0% (white). Labelling of wave vectors follows normal convention: T = [0,0,0], X = [1/2,0,0], M = [1/2,1/2,0] and 
R= [1/2,1/2,1/2]. 



C(t) = 



1 



N M-n 



N(M - n) 



Y^ Y v i ( mAt ) v J ( mAt + *) ( 3 ) 



where Vj is the velocity component of the jth atom. At 
is the time interval of 0.02 ps. M is the total number 
of time steps. The system was simulated for a total of 
50 ps which corresponds to M = 2500. The correlation 
function was calculated for each atom with a time length 
of 30 ps (n = 1500). A Gaussian profile was used be- 
fore Fourier transformation to suppress the ripple effect 
caused by time cut-off. The VACF of angular velocities of 
the C-N rigid rods rotating about their center of mass was 
also calculated. The Fourier transformation then gave us 
the DoS of the pure rotational modes of these rigid rods. 
To obtain the phonon spectrum, the calculated DoS was 
first multiplied by a weighting factor Airbk/mk contain- 
ing the scattering length b^ and the atomic mass m^ 
of the fcth atom. Then, in order to mimic experimental 
resolution^, the DoS is convolved with a Gaussian with 
FWHM of 10% of the energy transfer. Figure [7] shows the 
good agreement between the calculated spectrum and the 
experiment^ 

In order to draw links between the atomic vibrations 
and the phonon properties, we computed both the overall 
DoS and the angular DoS of rigid rod C-N, as shown in 
Figure |8ja) and (b), respectively. At zero pressure, the 
peak around 1.0 THz are stiffened with elevated tempera- 
ture as observed in the neutron-scattering experimenlpSl. 



We found that the temperature dependence of the peak 
energy is slightly sublinear. By comparing Figure p[b) 
with Figure (8ja), we found that part (about half) of the 
acoustic peak around 1.0 THz is from vibrations involv- 
ing rotations of the C-N rod around its centre of mass 
(bearing in mind that the angular DoS in Figure p[b) 
is renormalized, so that only the relative heights of the 
peaks are indicative). Optic peaks round 3.0 and 9.0 THz 
are hardly changed, because much of their motions can 
be cast onto the rotations of the C-N rod. The peak with 
the highest frequency in the overall DoS is from pure Zn- 
C/N bond flexing — it completely disappears in the an- 
gular DoS. Unlike the acoustic peak, this peak is softened 
on heating due to the thermal expansion of the Zn-C/N 
bond, and flattens with elevated temperature due to the 
finite life time of the corresponding phonons. 



We further calculated the DoS at ambient tempera- 
ture (300 K) under different pressures, as shown in Fig- 
ure [9] The acoustic peak around 1.0 THz is softened on 
compression. The optic peaks around 3.0 THz and 9.0 
THz follow the same trend, which means that all these 
modes have negative Griineisen parameters (7k, a) and 
contribute to the NTE of the material. The highest en- 
ergy peak around 14 THz is stiffened under compres- 
sion, suggesting positive 7k, a f° r the Zn-C/N bond flex- 
ing modes. 
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FIG. 11: (a) Total DoS from lattice-dynamics calculation 
shaded according to the average Griineisen parameter (7k. a) 
around each energy point. Modes with 7k,A less than —10 are 
in black. Modes with slightly negative 7k, a ~ —1 are in light 
grey. Positive-7k,A modes are in white, (b) Total DoS from 
lattice-dynamics calculation shaded according to how 'RUM' 
each mode is. Modes in black are RUMs and in white are 
non-RUMs. 



B. Rigid unit modes and negative thermal 
expansion 

To further investigate the nature of various peaks in 
the DoS, we decided to categorize the vibrational modes 
in the material using the rigid unit mode modeP^. First, 
we calculated 7k, a from phonon frequencies of expanded 
and contracted (±0.01%) unit-cell volumes. We then 
coloured the dispersion curves according to both mag- 
nitudes and signs of 7k, a j as shown in Figure 10 'a). This 



representation highlights the most important phonon 
branches responsible for NTE of the material, namely the 
low-lying acoustic modes around 1.0 THz and the lowest- 
energy optic branches around 2.0 THz, which both have 
the most negative 7k, a and span the entire Brillouin zone. 
This was also highlighted in the DoS in Figure [Tl|a), 
where we shaded the DoS according to the mean values 
°f 7k, a for each frequency bin. Then, we calculated the 
rigi d unit modes (RUMs) of Zn(CN) 2 using the CRUSH 
cod c 1 37 * 38 ^ as the set of eigenvectors of the dynamical ma- 
trix whose eigenvalues are zero. We took the dot product 
between the eigenvectors of the RUMs and the eigen- 
vectors from the lattice-dynamic calculation, and then 
coloured the dispersion curves in Figure [l0[ b) by the ex- 
tent to which each mode eigenvector can be described in 
terms of correlated whole-body translations (in blue) and 
rotations (in red) of [Zn(C/N)4] tetrahedra. The DoS in 
Figure fTlTb) was also shaded accordingly. One can see 
that all the modes with negative 7k, a that contribute to 
the NTE of the material are RUMs.' 



FIG. 12: Calculated Griineisen parameters (7k, a) of Zn(CN)2 
under different pressures (Black: 0.0 GPa; Blue: 0.2 GPa; Red: 
0.4 GPa). 7k,A of the modes around 0.5 THz, 2.0 THz and 
8.5 THz will become more negative with elevated pressure, 
resulting in a more negative 7 at higher pressure. 



we conclude that the acoustic modes around 1.0 THz, like 
those at M and X, are characterized by translational mo- 
tions of the rigid tetrahedral units, partly involving an- 
gular rotations of the C-N rod. The optic peaks around 
3.0 and 9.0 THz can be seen as the result of neighbouring 
tetrahedral units rotating against each other. The RUM 
nature of these modes guarantees their low frequencies 
and large negative 7k, a- The relatively high frequencies 
of the optic RUMs is due to the breaking of the Zn-C/N- 
N/C alignment and the magnitudes of their negative 7k,A 
suffer accordingly. Further study of the eigenvectors re- 
vealed that the non-RUM peaks between 9-10 THz cor- 
respond to the pure angular distortion of C/N-Zn-N/C 
within the tetrahedra. 

Values of the overall Griineisen parameter (7) calcu- 
lated using Eq. [2] are -1.15 at 0.0 GPa, -1.39 at 0.2 
GPa, -1.75 at 0.4 GPa and -2.36 at 0.6 GPa. Figure [12 
shows the profiles of 7k, a up to 0.4 GPa. It is clear that 7 
becomes more negative on compression due the contribu- 
tions from the RUMs around 1.0 THz, 3.0 THz and 9.0 
THz, with the translational RUMs having the most neg- 
ative 7k. a (refer to Figure H'a)) contributing the most. 
It turns out that, besides trie negative 7, Zn(CN) 2 also 
has negative d'y/dP and d 2; y/dP 2 . 

It is worth to mention at this point that the pressure 
dependence of 7 determines the pressure dependence of 
the coefficient of thermal expansion, ay, of the material. 
Follows from equation [T] one has 
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with B and B' the bulk modulus and its first deriva- 



Combining this with the findings in the former section, tive, respectively. Table III lists the values of the terms 



TABLE III: Computed values of the terms in equation [3] at different temperatures around P = using quasi- harmonic lattice 
dynamics (QHLD) and molecular dynamics (MD). Units of d In a/dP, (1 - B')/B, and (l/*f)(dy/dP) are GPa" 1 . Note that 
the values of B here, which obtained from fittings using Birch-Murnaghan equation of states based on the finite strain theory, 
are slightly different from the numerically calculated ones in Figure [5] 
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10 


-11.1 


-48.9 


1.77 


0.77 


45.15 


4.1 


-0.07 


1.44 


50 


-41.7 


-45.0 


0.88 


0.80 


43.51 


-3.0 


0.09 


0.98 


150 


-54.6 


-37.8 


0.78 


0.71 


39.81 


-6.0 


0.18 


0.93 


300 


-52.8 


-30.3 


0.85 


0.64 


34.47 


-4.2 


0.15 


1.01 


400 


-51.9 


-24.0 


0.87 


0.94 


32.22 


-3.1 


0.13 


1.03 


500 


-51.0 


-18.9 


0.89 


1.11 


30.11 


-1.8 


0.09 


1.05 



"Calculated in QHLD. 

'Obtained using the 3rd-order Birch-Murnaghan equation of 
states to fit to the isotherms from the MD. 



in this equation at different temperatures. It is possi- 
ble to use both MD and QHLD to calculate values of 
d In ay/dP — differences are to be expected because of 
anharmonic effects, and due to the fact that the QHLD 
method uses quantum statistics but the MD is classical 
- it is only possible to calculate 7 and dlnj/dP from 
QHLD. Because of these differences, it is unlikely that we 
can obtain perfect agreement between both sides of the 
equation. However, the data in the table clearly indicate 
that the change of 7 with pressure contribute dominantly 
to the pressure dependence of ay 



VIII. STRUCTURAL STUDY 

A. The Local Picture on Compression and Heating 

In the real space, variations of geometrical features 
of the material, such as the Zn-C/N bond length, the 
N/C-Zn-C/N angle within the tetrahedral unit, and the 
Zn-C/N-N/C angle are important for us to build a local 
picture of the system subject to both compression and 
heating. 

Distributions of these quantities, as shown in Fig- 
were obtained from the atomic trajectory data 
MD. The large spread of the distributions sug- 
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of the 

gests large vibrations of these quantities at high tem- 
perature. The slightly expansion-biased broadening of 
the bond-length distribution on heating indicates an en- 
hanced thermal expansion in the bond. We found that 
the deviation of the average N/C-Zn-C/N angle in the 
tetrahedral unit from its equilibrium of 109.47° is triv- 
ially small (~ 0.2°) even at very high temperature (~ 600 
K) and pressure (~ 0.6 GPa). 

The average bond length and the average angle dis- 
tortion in Zn(CN)2 as functions of both temperature 
and pressure are shown in Figure |14| Compression pro- 
gressively increases the Zn-C/N-N/C angle with elevated 
temperature. At very low and zero temperature, the 



trend is that the angle will be hardly changed by pres- 
sure. This is exactly what we expected in the previous 
section where we suggested that, at zero temperature, 
the volume change of the material due to the pressure 
arises solely from the compression of the Zn-C/N bonds, 
resulting in positive B' . It is when the Zn-C/N-N/C an- 
gle starts to increase under compression and contribute 
to the volume change that the material shows negative 
B' , i.e. the pressure- induced softening. The inset of Fig- 
ure 14 '&) shows the increase of bond length with temper- 



ature. The superlinear behaviour indicates the softening 
of the bond at higher temperature due to the thermal ex- 
pansion. In Figure |l4[ b) , the inset shows the increase of 
the Zn-C/N-N/C angle with elevated temperature. The 
sublinear behaviour suggests that the angle will become 
more rigid, resulting in a less negative B' at high tem- 
peratures. Both the superlinearity and the sublinearity 
in the plots rely on the capture of the anharmonicity of 
the material. 

To see the real-space picture of RUMs in Zn(CN)2, 
we quantified the proportion of thermal-excited rigid- 
unit rotations at different pressures. Using our GASP 
code^2M2 we compared ten snapshots of an MD simu- 
lation, where each snapshot is separated by 2 ps, with 
the ideal structure. This is repeated for various pressures 
and temperatures. GASP, using geometric algebra, can 
partition the atomic displacements for every comparison 
made into the mean squared rigid-tetrahedron rotations, 
translational displacements and unit deformations. We 
then computed the average proportion of rigid-unit ro- 
tations at each temperature. To exclude those rigid-unit 
rotations due to pure topology reasons, i.e. rotations that 
accidentally maintain the shape of the tetrahedral unit 
under thermal excitation but are not because of the fea- 
tures of motion, we set up a benchmark calculation using 
an ideal cristobalite structure without any interactions 
other than bonds to hold Zn-C/N and C-N. The reason 
to use the single- framework lattice is to avoid the problem 
of two interpenetrating frameworks crushing into each 
other in the MD due to the lack of long-distance interac- 
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FIG. 13: Distributions of (a) the Zn-C/N bond length; (b) the 
cosine of Zn-C/N-N/C angle distortion; (c) the N/C-Zn-C/N 
angle in the tetrahedral unit at 0.0 GPa. 



tions. Figure [15] shows the results at different pressures 
and temperatures with coloured areas. The proportion of 
the rigid-unit rotation of a real silica system is given in 
the plot as a comparison. The figure suggests that com- 
pression will enhance the rigid-unit rotations. At ambi- 
ent temperature ~ 300 K, for example, the average pro- 
portions are 63% at 0.0 GPa, 64% at 0.2 GPa, 65% at 
0.4 GPa and 66% at 0.6 GPa, compared to 35% of the 
benchmark and 90% of the silica system. Another im- 
portant point is that, at certain pressure, the proportion 
of rigid-unit rotation will decrease with elevated temper- 
ature. This trend corresponds to the peaks around 1.0, 
3.0 and 9.0 THz in DoS stiffened on heating, as shown 
in Figure [8] The Griineisen parameters of these modes 
will consequently become less negative, and so does the 
coefficient of thermal expansion. 



B. The Nearest Neighbouring Zn. . .Zn Distance 

There is a puzzling observation from the X-ray pair 
distribution function measurements of Zn(CN)^that in- 
stantaneous Zn. . .Zn distances contract less rapidly on 
heating than does the cell length, while crystallographi- 
cally the two should be linked. The explanation relies on 
the nature of the acoustic modes as translational RUMs. 

As mentioned before, the acoustic modes around 1.0 
THz are translational RUMs. With the absence of these 
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FIG. 14: (a) The average Zn-C/N bond length increases with 
elevated temperature. The inset shows superlinear behaviour 
of the bond length on heating, (b) The average Zn-C/N-N/C 
angle distortion increases with elevated temperature and pres- 
sure. The inset shows the sublinear behavior of the angle dis- 
tortion on heating. 



modes, one can estimate 7 to be —0.53 compared to 
— 1.15 when all modes are included. According to Eq. [T] 
this tells that the translational RUMs count for about 
half the NTE of the material. Unlike rotational RUMs at 
3.0 and 9.0 THz that will reduce the nearest Zn. . .Zn dis- 
tance, the translational RUMs correspond to collective 
translations of the neighbouring rigid units that moves 
zinc atoms off site and retains the distance of the nearest 
neighbouring zincs. This kind of vibration involves the 
rotation of C-N rod around its centre of mass, consistent 
with the previous finding in the DoS that part of the 
acoustic modes are from the C-N rod rotations. 

We calculated the average distance of the nearest 
neighbouring Zn. . .Zn using the the trajectory data from 
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FIG. 15: Computed proportion of rigid-unit rotations under 
thermal excitation by using GASP to analyze the MD data. 
Boundary of areas with different colors represents the value 
of proportion at different pressure: yellow for 0.0 GPa, orange 
for 0.2 GPa, pink for 0.4 GPa and red for 0.6 GPa. Green 
area and dark-red area are for the benchmark calculation and 
silica system, respectively. The results clearly show the en- 
hancement of rigid-unit rotations hence the rotational RUMs 
under pressure. 



the MD simulations. Figure [16] shows the temperature 
dependence of the average distance under different pres- 
sures. At zero pressure, the ratio between the linear coeffi- 
cient of thermal expansion (CTE) of the average Zn. . .Zn 
distance, ctznZm an d the overall linear CTE of the ma- 
terial, a, is 0.67, with aznZn = — 8.5MK -1 and a — 
— 12.62 MK" 1 . This result agrees well with the experi- 
mental raticP of 0.71, with a Zn Zn = -12.42(12) MK" 1 
and a = -17.40(18) MK' 1 . 



IX. DISCUSSIONS AND CONCLUSIONS 

Although the increase of pressure or temperature 
would both result in volume contraction in Zn(CN)2, 
the pressure and temperature dependence of the NTE 
in Zn(CN)2 are totally different. The two variables have 
the opposite effects on the normal modes of the material. 
Increasing temperature stiffens the low-frequency peaks 
and softens the high-frequency peaks in the DoS (Fig- 
urepj) accompanied by the reduction of NTE, while com- 
pression would soften the low-frequency peaks and stiffen 
the high-frequency peaks [Figure [9] resulting in NTE en- 
hancement. Raising temperature slows the mode soft- 
ening caused by compression, and postpones the phase 
transition. The enormous decrease of the bulk modulus 
on heating contrasts the small change of the bulk modu- 
lus on compression, which disobeys Birch's law of corre- 
sponding states. 

The pressure and temperature dependence of the Zn- 



FIG. 16: Calculated average distances for the nearest- 
neighboring Zn. . .Zn as functions of temperature at different 
pressures. The inset shows the coefficients of thermal expan- 
sion of the distance at different pressures averaged over the 
temperature range. 



C/N bond and the N/C-Zn-C/N angle are also intriguing. 
The large vibrational amplitude of the N/C-Zn-C/N an- 
gle extends the Zn-C/N bond, which can be seen (under 
constant pressure) in Figure 
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The superlinearity sug- 
gests an enhancement of the thermal expansion in the 
bond with more distorted angle at higher temperature. 
The sublinearity of the average angle distortion corre- 
sponds to the stiffened RUMs involving the rotational 
vibrations of the C-N rods. 

The ability to carry out MD for Zn(CN)2 using the 
potential model is vital in this study. It allows us to cap- 
ture all of the anharmonicity, to reproduce all of the ex- 
otic properties of the material, and to study their pres- 
sure and temperature dependence. The origins of various 
properties are revealed by linking features in both energy 
and real space. 

One significant example is the pressure-enhanced NTE, 
which has been well reproduced in Figure [3|a) . This be- 
haviour is linked to a feature in energy space that the 
peaks around 1.0 and 3.0 THz are softened under com- 
pression (as in Figure p3J), as well as the rising proportion 
of the rigid-unit rotations in real space on compression 
(as in Figure 15). 

Another example is the reduction of NTE with elevated 
temperature. Stiffening of the modes around 1.0 and 3.0 
THz on heating (as shown in the DoS in Figure pi) makes 
their Griineisen parameters less negative, hence the re- 
duction of NTE in the material. In real space, Figure [TB] 
clearly shows the trend of decrease in the proportion of 
rigid-unit rotations on heating. 

The third example is the temperature dependence of 
the bulk modulus and its first derivative B' . The large 
decrease of the bulk modulus on heating is due to the 
involvement of the Zn-C/N-N/C angular vibrations at 
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non-zero temperature and the softening of the Zn-C/N 
bond, which can be seen both in the DoS (Figure^ a)) 
and in the real-space picture of Figure [13] and Figure [14] 
As shown in Figure [l4|b), the average Zn-C/N-N/C an- 
gle starts to possess an almost linear pressure dependence 
at medium temperature and contributes to the pressure- 
induced volume change of the material. According to 
a simple geometrical relation^, the relative decrease in 
cell is roughly proportional to the Zn-C/N-N/C angle 
squared, hence to the pressure squared, resulting in neg- 
ative Bq. However, this contribution from the pressure- 
induced change in the Zn-C/N-N/C angle to the volume 
contraction will be hindered at high temperature due to 
the anharmonic stiffening of the corresponding modes 
(Figure [81a)). On the other hand, the contribution from 
Zn-C/N bond compression will increase due to the ther- 
mal softening of the bond. The combined effect of these 
two trends is expected to result in a less negative B' at 
high temperatures (see values in Table III ) . 
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